﻿///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// TetGen                                                                    //
//                                                                           //
// A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator       //
//                                                                           //
// Version 1.5                                                               //
// November 4, 2013                                                          //
//                                                                           //
// TetGen is freely available through the website: http://www.tetgen.org.    //
//   It may be copied, modified, and redistributed for non-commercial use.   //
//   Please consult the file LICENSE for the detailed copyright notices.     //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////

#pragma once
// TetGen only uses the C standard library.
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <ctime>
#include <cassert> 
#include <iostream>
#include <string>

using std::cout;
using std::endl;
using std::string;

namespace Tetgen
{
	/*
	 To compile TetGen as a library instead of an executable program, define
	   the TETLIBRARY symbol.

	 #define TETLIBRARY

	 Uncomment the following line to disable assert macros. These macros were
	   inserted in the code where I hoped to catch bugs. They may slow down the
	   speed of TetGen.

	 #define NDEBUG

	 TetGen default uses the double precision (64 bit) for a real number. 
	   Alternatively, one can use the single precision (32 bit) 'float' if the
	   memory is limited.
	*/

	typedef double REAL;  //typedef float REAL

	// Maximum number of characters in a file name (including the null).文件名最大长度
	const int FILENAMESIZE = 1024;

	// Maximum number of chars in a line read from a file (including the null).一次从文件中读取的最大长度
	const int INPUTLINESIZE = 2048;

	/* The types 'intptr_t' and 'uintptr_t' are signed and unsigned integer types,
	//   respectively. They are guaranteed to be the same width as a pointer.
	//   They are defined in <stdint.h> by the C99 Standard. However, Microsoft 
	//   Visual C++ 2003 -- 2008 (Visual C++ 7.1 - 9) doesn't ship with this header
	//   file. In such case, we can define them by ourself. 
	// Update (learned from Stack Overflow): Visual Studio 2010 and Visual C++ 2010
	//   Express both have stdint.h

	// The following piece of code was provided by Steven Johnson (MIT). Define the
	//   symbol _MSC_VER if you are using Microsoft Visual C++. Moreover, define 
	//   the _WIN64 symbol if you are running TetGen on Win64 systems.*/

#ifdef _MSC_VER // Microsoft Visual C++
#  ifdef _WIN64
	typedef __int64 intptr_t;
	typedef unsigned __int64 uintptr_t;
#  else // not _WIN64
	typedef int intptr_t;
	typedef unsigned int uintptr_t;
#  endif
#else // not Visual C++
#  include <stdint.h>
#endif

	// The input object of TetGen. They are recognized by either the input 
	//   file extensions or by the specified options. 
	// Currently the following objects are supported:
	//   - NODES, a list of nodes (.node); 
	//   - POLY, a piecewise linear complex (.poly or .smesh); 
	//   - OFF, a polyhedron (.off, Geomview's file format); 
	//   - PLY, a polyhedron (.ply, file format from gatech, only ASCII);
	//   - STL, a surface mesh (.stl, stereolithography format);
	//   - MEDIT, a surface mesh (.mesh, Medit's file format); 
	//   - MESH, a tetrahedral mesh (.ele).
	// If no extension is available, the imposed command line switch
	//   (-p or -r) implies the object. 
	enum class TetInEnum
	{
		NODES,
		POLY,
		OFF,
		PLY,
		STL,
		MEDIT,
		VTK,
		MESH
	};


	/*//////////////////////////////////////////////////////////////////////////////
	//                                                                           //
	// tetgenio                                                                  //
	// 为转换网格定义的数据结构                                                                          //
	// A structure for transferring data into and out of TetGen's mesh structure,//
	// 'tetgenmesh' (declared below).                                            //
	//                                                                           //
	// The input of TetGen is either a 3D point set, or a 3D piecewise linear    //
	// complex (PLC), or a tetrahedral mesh.  Depending on the input object and  //
	// the specified options, the output of TetGen is either a Delaunay (or wei- //
	// ghted Delaunay) tetrahedralization, or a constrained (Delaunay) tetrahed- //
	// ralization, or a quality tetrahedral mesh.                                //
	//                                                                           //
	// A piecewise linear complex (PLC) represents a 3D polyhedral domain with   //
	// possibly internal boundaries(subdomains). It is introduced in [Miller et  //
	// al, 1996]. Basically it is a set of "cells", i.e., vertices, edges, poly- //
	// gons, and polyhedra, and the intersection of any two of its cells is the  //
	// union of other cells of it.                                               //
	//                                                                           //
	// TetGen uses a set of files to describe the inputs and outputs. Each file  //
	// is identified from its file extension (.node, .ele, .face, .edge, etc).   //
	//                                                                           //
	// The 'tetgenio' structure is a collection of arrays of data, i.e., points, //
	// facets, tetrahedra, and so forth. It contains functions to read and write //
	// (input and output) files of TetGen as well as other supported mesh files. //
	//                                                                           //
	// Once an object of tetgenio is declared,  no array is created. One has to  //
	// allocate enough memory for them. On deletion of this object, the memory   //
	// occupied by these arrays needs to be freed.  The routine deinitialize()   //
	// will be automatically called.  It frees the memory for an array if it is  //
	// not a NULL. Note that it assumes that the memory is allocated by the C++  //
	// "new" operator.  Otherwise, the user is responsible to free them and all  //
	// pointers must be NULL before the call of the destructor.                  //
	//                                                                           //
	//////////////////////////////////////////////////////////////////////////////*/
	class tetgenio
	{
	public:

		// A "polygon" describes a simple polygon (no holes). It is not necessarily
		//   convex. Each polygon contains a number of corners (points) and the same
		//   number of sides (edges).  The points of the polygon must be given in
		//   either counterclockwise or clockwise order and they form a ring, so 
		//   every two consecutive points forms an edge of the polygon.
		struct  polygon		//面片
		{
			int *vertexlist;
			int numberofvertices;
		};

		// A "facet" describes a polygonal region possibly with holes, edges, and 
		//   points floating in it.  Each facet consists of a list of polygons and
		//   a list of hole points (which lie strictly inside holes).面元
		struct facet
		{
			polygon *polygonlist;
			int numberofpolygons;
			REAL *holelist;
			int numberofholes;
		};

		// A "voroedge" is an edge of the Voronoi diagram. It corresponds to a
		//   Delaunay face.  Each voroedge is either a line segment connecting
		//   two Voronoi vertices or a ray starting from a Voronoi vertex to an
		//   "infinite vertex".  'v1' and 'v2' are two indices pointing to the
		//   list of Voronoi vertices. 'v1' must be non-negative, while 'v2' may
		//   be -1 if it is a ray, in this case, the unit normal of this ray is
		//   given in 'vnormal'. 
		struct voroedge
		{
			int v1, v2;
			REAL vnormal[3];
		};

		// A "vorofacet" is an facet of the Voronoi diagram. It corresponds to a
		//   Delaunay edge.  Each Voronoi facet is a convex polygon formed by a
		//   list of Voronoi edges, it may not be closed.  'c1' and 'c2' are two
		//   indices pointing into the list of Voronoi cells, i.e., the two cells
		//   share this facet.  'elist' is an array of indices pointing into the
		//   list of Voronoi edges, 'elist[0]' saves the number of Voronoi edges
		//   (including rays) of this facet.
		struct vorofacet
		{
			int c1, c2;
			int *elist;
		} ;


		// Additional parameters associated with an input (or mesh) vertex.
		//   These informations are provided by CAD libraries. 
		struct pointparam
		{
			REAL uv[2];
			int tag;
			int type; // 0, 1, or 2.
		};

		// Callback functions for meshing PSCs.
		typedef REAL(*GetVertexParamOnEdge)(void*, int, int);
		typedef void(*GetSteinerOnEdge)(void*, int, REAL, REAL*);
		typedef void(*GetVertexParamOnFace)(void*, int, int, REAL*);
		typedef void(*GetEdgeSteinerParamOnFace)(void*, int, REAL, int, REAL*);
		typedef void(*GetSteinerOnFace)(void*, int, REAL*, REAL*);

		// A callback function for mesh refinement.
		typedef bool(*TetSizeFunc)(REAL*, REAL*, REAL*, REAL*, REAL*, REAL);

		// Items are numbered starting from 'firstnumber' (0 or 1), default is 0.起始编号，默认从0开始
		int firstnumber=0;

		// Dimension of the mesh (2 or 3), default is 3.网格维度
		int mesh_dim=3;

		// Does the lines in .node file contain index or not, default is true.节点文件是否包含编号
		bool useindex=true;

		/*
		 'pointlist':  An array of point coordinates.  The first point's x
		   coordinate is at index [0] and its y coordinate at index [1], its
		   z coordinate is at index [2], followed by the coordinates of the
		   remaining points.  Each point occupies three REALs. 
		 'pointattributelist':  An array of point attributes.  Each point's
		   attributes occupy 'numberofpointattributes' REALs.
		 'pointmtrlist': An array of metric tensors at points. Each point's
		   tensor occupies 'numberofpointmtr' REALs.
		 'pointmarkerlist':  An array of point markers; one integer per point.
		*/
		REAL *pointlist;
		REAL *pointattributelist;
		REAL *pointmtrlist;
		int  *pointmarkerlist;
		pointparam *pointparamlist;
		int numberofpoints=0;
		int numberofpointattributes=0;
		int numberofpointmtrs;

		/*
		 'tetrahedronlist':  An array of tetrahedron corners.  The first 
		   tetrahedron's first corner is at index [0], followed by its other 
		   corners, followed by six nodes on the edges of the tetrahedron if the
		   second order option (-o2) is applied. Each tetrahedron occupies
		   'numberofcorners' ints.  The second order nodes are ouput only. 
		 'tetrahedronattributelist':  An array of tetrahedron attributes.  Each
		   tetrahedron's attributes occupy 'numberoftetrahedronattributes' REALs.
		 'tetrahedronvolumelist':  An array of constraints, i.e. tetrahedron's
		   volume; one REAL per element.  Input only.
		 'neighborlist':  An array of tetrahedron neighbors; 4 ints per element. 
		   Output only.
		*/
		int  *tetrahedronlist;
		REAL *tetrahedronattributelist;
		REAL *tetrahedronvolumelist;
		int  *neighborlist;
		int numberoftetrahedra;
		int numberofcorners;
		int numberoftetrahedronattributes;

		// 'facetlist':  An array of facets.  Each entry is a structure of facet.
		// 'facetmarkerlist':  An array of facet markers; one int per facet.
		facet *facetlist;
		int *facetmarkerlist;
		int numberoffacets;

		// 'holelist':  An array of holes (in volume).  Each hole is given by a
		//   seed (point) which lies strictly inside it. The first seed's x, y and z
		//   coordinates are at indices [0], [1] and [2], followed by the
		//   remaining seeds.  Three REALs per hole. 
		REAL *holelist;
		int numberofholes;

		// 'regionlist': An array of regions (subdomains).  Each region is given by
		//   a seed (point) which lies strictly inside it. The first seed's x, y and
		//   z coordinates are at indices [0], [1] and [2], followed by the regional
		//   attribute at index [3], followed by the maximum volume at index [4]. 
		//   Five REALs per region.
		// Note that each regional attribute is used only if you select the 'A'
		//   switch, and each volume constraint is used only if you select the
		//   'a' switch (with no number following).
		REAL *regionlist;
		int numberofregions;

		// 'facetconstraintlist':  An array of facet constraints.  Each constraint
		//   specifies a maximum area bound on the subfaces of that facet.  The
		//   first facet constraint is given by a facet marker at index [0] and its
		//   maximum area bound at index [1], followed by the remaining facet con-
		//   straints. Two REALs per facet constraint.  Note: the facet marker is
		//   actually an integer.
		REAL *facetconstraintlist;
		int numberoffacetconstraints;

		// 'segmentconstraintlist': An array of segment constraints. Each constraint 
		//   specifies a maximum length bound on the subsegments of that segment.
		//   The first constraint is given by the two endpoints of the segment at
		//   index [0] and [1], and the maximum length bound at index [2], followed
		//   by the remaining segment constraints.  Three REALs per constraint. 
		//   Note the segment endpoints are actually integers.
		REAL *segmentconstraintlist;
		int numberofsegmentconstraints;

		// 'trifacelist':  An array of face (triangle) corners.  The first face's
		//   three corners are at indices [0], [1] and [2], followed by the remaining
		//   faces.  Three ints per face.
		// 'trifacemarkerlist':  An array of face markers; one int per face.
		// 'o2facelist':  An array of second order nodes (on the edges) of the face.
		//   It is output only if the second order option (-o2) is applied. The
		//   first face's three second order nodes are at [0], [1], and [2],
		//   followed by the remaining faces.  Three ints per face.
		// 'adjtetlist':  An array of adjacent tetrahedra to the faces. The first
		//   face's two adjacent tetrahedra are at indices [0] and [1], followed by
		//   the remaining faces.  A '-1' indicates outside (no adj. tet). This list
		//   is output when '-nn' switch is used. Output only.
		int *trifacelist;
		int *trifacemarkerlist;
		int *o2facelist;
		int *adjtetlist;
		int numberoftrifaces;

		// 'edgelist':  An array of edge endpoints.  The first edge's endpoints
		//   are at indices [0] and [1], followed by the remaining edges.
		//   Two ints per edge.
		// 'edgemarkerlist':  An array of edge markers; one int per edge.
		// 'o2edgelist':  An array of midpoints of edges. It is output only if the
		//   second order option (-o2) is applied. One int per edge.
		// 'edgeadjtetlist':  An array of adjacent tetrahedra to the edges.  One
		//   tetrahedron (an integer) per edge.
		int *edgelist;
		int *edgemarkerlist;
		int *o2edgelist;
		int *edgeadjtetlist;
		int numberofedges;

		// 'vpointlist':  An array of Voronoi vertex coordinates (like pointlist).
		// 'vedgelist':  An array of Voronoi edges.  Each entry is a 'voroedge'.
		// 'vfacetlist':  An array of Voronoi facets. Each entry is a 'vorofacet'.
		// 'vcelllist':  An array of Voronoi cells.  Each entry is an array of
		//   indices pointing into 'vfacetlist'. The 0th entry is used to store
		//   the length of this array.
		REAL *vpointlist;
		voroedge *vedgelist;
		vorofacet *vfacetlist;
		int **vcelllist;
		int numberofvpoints;
		int numberofvedges;
		int numberofvfacets;
		int numberofvcells;

		// Variable (and callback functions) for meshing PSCs.
		void *geomhandle;
		GetVertexParamOnEdge getvertexparamonedge;
		GetSteinerOnEdge getsteineronedge;
		GetVertexParamOnFace getvertexparamonface;
		GetEdgeSteinerParamOnFace getedgesteinerparamonface;
		GetSteinerOnFace getsteineronface;

	public:
		static void init(polygon* p);
		static void init(facet* f);

		// Constructor & destructor.
		tetgenio();
		~tetgenio();

		// A callback function.
		TetSizeFunc tetunsuitable;

		// Input & output routines.
		bool load_node_call(FILE* infile, int markers, int uvflag, char*);
		bool load_node(const string& fname);
		bool load_edge(char*);
		bool load_face(char*);
		bool load_tet(char*);
		bool load_vol(char*);
		bool load_var(char*);
		bool load_mtr(char*);
		bool load_poly(const string& fname);
		bool load_off(char*);
		bool load_ply(char*);
		bool load_stl(char*);
		bool load_vtk(char*);
		bool load_medit(char*, int);
		bool load_plc(char*, const TetInEnum& ob);
		bool load_tetmesh(char*, const TetInEnum& object);
		void save_nodes(char*);
		void save_elements(char*);
		void save_faces(char*);
		void save_edges(char*);
		void save_neighbors(char*);
		void save_poly(char*);
		void save_faces2smesh(char*);

		// Read line and parse string functions.
		char *readline(char* string, FILE* infile, int *linenumber);
		char *findnextfield(char* string);
		char *readnumberline(char* string, FILE* infile, char* infilename);
		char *findnextnumber(char* string);
	};

	/*//////////////////////////////////////////////////////////////////////////////
	//	网格控制参数                                                                           //
	// tetgenbehavior                                                            //
	//                                                                           //
	// A structure for maintaining the switches and parameters used by TetGen's  //
	// mesh data structure and algorithms.                                       //
	//                                                                           //
	// All switches and parameters are initialized with default values. They can //
	// be set by the command line arguments (a list of strings) of TetGen.       //
	//                                                                           //
	// NOTE: Some of the switches are incompatible. While some may depend on     //
	// other switches.  The routine parse_commandline() sets the switches from   //
	// the command line (a list of strings) and checks the consistency of the    //
	// applied switches.                                                         //
	//                                                                           //
	//////////////////////////////////////////////////////////////////////////////*/
	class tetgenbehavior
	{
	public:
		tetgenbehavior();
		virtual ~tetgenbehavior();

		void syntax();
		void usage();

		void setDefaultOption();	//gzz默认选项
		void copyFileName(const tetgenbehavior& b);

		// Command line parse routine.
		bool parse_commandline(int argc, char **argv);
		bool parse_commandline(char *switches);

		bool Refine() const;
		bool Quiet() const;
		bool Plc() const;
		void setPlc(bool val);

		REAL MinDihedral() const;
		REAL MaxVolume() const;
		REAL MinRatio() const;

		void setMinDihedral(double val);
		void setMaxVolume(double val);
		void setMinRatio(double val);

	public:
		// Switches of TetGen. 
		int psc=0;                                                         // '-s', 0.
		int quality;                                                     // '-q', 0.
		int nobisect;                                                    // '-Y', 0.
		int coarsen;                                                     // '-R', 0.
		int weighted;                                                    // '-w', 0.
		int brio_hilbert;                                                // '-b', 1.
		int incrflip;                                                    // '-l', 0.
		int flipinsert;                                                  // '-L', 0.
		int metric;                                                      // '-m', 0.
		int varvolume;                                                   // '-a', 0.
		int fixedvolume;                                                 // '-a', 0.
		int regionattrib;                                                // '-A', 0.
		int conforming;                                                  // '-D', 0.
		int insertaddpoints;                                             // '-i', 0.
		int diagnose;                                                    // '-d', 0.
		int convex;                                                      // '-c', 0.
		int nomergefacet;                                                // '-M', 0.
		int nomergevertex;                                               // '-M', 0.
		int noexact;                                                     // '-X', 0.
		int nostaticfilter;                                              // '-X', 0.
		int zeroindex;                                                   // '-z', 0.
		int facesout;                                                    // '-f', 0.
		int edgesout;                                                    // '-e', 0.
		int neighout;                                                    // '-n', 0.
		int voroout;                                                     // '-v', 0.
		int meditview;                                                   // '-g', 0.
		int vtkview;                                                     // '-k', 0.

		int nobound;                                                     // '-B', 0.
		int nonodewritten;                                               // '-N', 0.
		int noelewritten;                                                // '-E', 0.
		int nofacewritten;                                               // '-F', 0.

		int noiterationnum;                                              // '-I', 0.
		int nojettison;                                                  // '-J', 0.
		int reversetetori;                                               // '-R', 0.
		int docheck;                                                     // '-C', 0.
		int verbose;                                                     // '-V', 0.

		// Parameters of TetGen. 
		int vertexperblock;                                           // '-x', 4092.
		int tetrahedraperblock;                                       // '-x', 8188.
		int shellfaceperblock;                                        // '-x', 2044.
		int nobisect_param;                                              // '-Y', 2.
		int addsteiner_algo;                                            // '-Y/', 1.
		int coarsen_param;                                               // '-R', 0.
		int weighted_param;                                              // '-w', 0.
		int fliplinklevel;                                                    // -1.
		int flipstarsize;                                                     // -1.
		int fliplinklevelinc;                                                 //  1.
		int reflevel;                                                    // '-D', 3.
		int optlevel;                                                    // '-O', 2.
		int optscheme;                                                   // '-O', 7.
		int delmaxfliplevel;                                                   // 1.
		int order;                                                       // '-o', 1.
		int steinerleft;                                                 // '-S', 0.
		int no_sort;                                                           // 0.
		int hilbert_order;                                           // '-b///', 52.
		int hilbert_limit;                                             // '-b//'  8.
		int brio_threshold;                                              // '-b' 64.

		REAL brio_ratio;                                             // '-b/' 0.125.
		REAL facet_ang_tol;                                          // '-p', 179.9.
		REAL optmaxdihedral;                                               // 165.0.
		REAL optminsmtdihed;                                               // 179.0.
		REAL optminslidihed;                                               // 179.0.  
		REAL epsilon;                                               // '-T', 1.0e-8.
		REAL minedgelength;                                                  // 0.0.
		REAL coarsen_percent;                                         // -R1/#, 1.0.

		// Strings of command line arguments and input/output file names.
		char commandline[1024];
		char infilename[1024];
		char outfilename[1024];
		char addinfilename[1024];
		char bgmeshfilename[1024];

		//input file type enum
		TetInEnum object = TetInEnum::POLY;

		//输出文件名称
		//string outfilename;
	protected:
		bool plc = false;                                                         // '-p', 0.
		bool refine = false;                                                      // '-r', 0.
		bool quiet = false;                                                       // '-Q', 0.

		//四面体单元的最小二面角
		REAL minDihedral=5.0;                                              // '-q', 5.0.
		//单元的最大体积限制
		REAL maxvolume=-1.0;                                               // '-a', -1.0.
		//单元的径边比
		REAL minratio=0.0;                                                 // '-q', 0.0.
	};

	///////////////////////////////////////////////////////////////////////////////
	// tetgenmesh                                                                //
	//                                                                           //
	// A structure for creating and updating tetrahedral meshes.                 //
	// 网格数据结构                                                                        //
	///////////////////////////////////////////////////////////////////////////////
	class tetgenmesh
	{
	public:
		tetgenmesh();
		~tetgenmesh();

		void freememory();

		void tetrahedralize(tetgenbehavior& b, tetgenio *in, tetgenio *out, tetgenio *addin = NULL, tetgenio *bgmin = NULL);
		void putGuoMesh(const int& cpart);

	public:

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh data structure                                                       //
		//                                                                           //
		// A tetrahedral mesh T of a 3D piecewise linear complex (PLC) X is a 3D     //
		// simplicial complex whose underlying space is equal to the space of X.  T  //
		// contains a 2D subcomplex S which is a triangular mesh of the boundary of  //
		// X. S contains a 1D subcomplex L which is a linear mesh of the boundary of //
		// S. Faces and edges in S and L are respectively called subfaces and segme- //
		// nts to distinguish them from others in T.                                 //
		//                                                                           //
		//TetGen stores the tetrahedra and vertices of T. 
		//The basic structure of a tetrahedron contains pointers to its vertices and adjacent tetrahedra. 
		//A vertex stores its x-, y-, and z-coordinates, and a pointer to a tetrahedron containing it. 
		//Both tetrahedra and vertices may contain user data.
		//                                                                           //
		// Each face of T belongs to either two tetrahedra or one tetrahedron. In    //
		// the latter case, the face is an exterior boundary face of T.  TetGen adds //
		// fictitious tetrahedra (one-to-one) at such faces, and connects them to an //
		// "infinite vertex" (which has no geometric coordinates).  One can imagine  //
		// such a vertex lies in 4D space and is visible by all exterior boundary    //
		// faces.  The extended set of tetrahedra (including the infinite vertex) is //
		// a tetrahedralization of a 3-pseudomanifold without boundary.  It has the  //
		// property that every face is shared by exactly two tetrahedra.             // 
		//                                                                           //
		// The current version of TetGen stores explicitly the subfaces and segments //
		// (which are in surface mesh S and the linear mesh L), respectively.  Extra //
		// pointers are allocated in tetrahedra and subfaces to point each others.   //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// The tetrahedron data structure.  It includes the following fields:
		//   - a list of four adjoining tetrahedra;
		//   - a list of four vertices;
		//   - a pointer to a list of four subfaces (optional, for -p switch);
		//   - a pointer to a list of six segments  (optional, for -p switch);
		//   - a list of user-defined floating-point attributes (optional);
		//   - a volume constraint (optional, for -a switch);
		//   - an integer of element marker (and flags);
		// The structure of a tetrahedron is an array of pointers.  Its actual size
		//   (the length of the array) is determined at runtime.

		typedef REAL **tetrahedron;

		// The subface data structure.  It includes the following fields:
		//   - a list of three adjoining subfaces;
		//   - a list of three vertices;
		//   - a list of three adjoining segments;
		//   - two adjoining tetrahedra;
		//   - an area constraint (optional, for -q switch);
		//   - an integer for boundary marker;
		//   - an integer for type, flags, etc.

		typedef REAL **shellface;

		// The point data structure.  It includes the following fields:
		//   - x, y and z coordinates;
		//   - a list of user-defined point attributes (optional);
		//   - u, v coordinates (optional, for -s switch);
		//   - a metric tensor (optional, for -q or -m switch);
		//   - a pointer to an adjacent tetrahedron;
		//   - a pointer to a parent (or a duplicate) point;
		//   - a pointer to an adjacent subface or segment (optional, -p switch);
		//   - a pointer to a tet in background mesh (optional, for -m switch);
		//   - an integer for boundary marker (point index);
		//   - an integer for point type (and flags).
		//   - an integer for geometry tag (optional, for -s switch).
		// The structure of a point is an array of REALs.  Its acutal size is 
		//   determined at the runtime.

		typedef REAL *point;

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Handles                                                                   //
		//                                                                           //
		// Navigation and manipulation in a tetrahedralization are accomplished by   //
		// operating on structures referred as ``handles". A handle is a pair (t,v), //
		// where t is a pointer to a tetrahedron, and v is a 4-bit integer, in the   //
		// range from 0 to 11. v is called the ``version'' of a tetrahedron, it rep- //
		// resents a directed edge of a specific face of the tetrahedron.            //
		//                                                                           //
		// There are 12 even permutations of the four vertices, each of them corres- //
		// ponds to a directed edge (a version) of the tetrahedron.  The 12 versions //
		// can be grouped into 4 distinct ``edge rings'' in 4 ``oriented faces'' of  //
		// this tetrahedron.  One can encode each version (a directed edge) into a   //
		// 4-bit integer such that the two upper bits encode the index (from 0 to 2) //
		// of this edge in the edge ring, and the two lower bits encode the index (  //
		// from 0 to 3) of the oriented face which contains this edge.               //  
		//                                                                           //
		// The four vertices of a tetrahedron are indexed from 0 to 3 (according to  //
		// their storage in the data structure).  Give each face the same index as   //
		// the node opposite it in the tetrahedron.  Denote the edge connecting face //
		// i to face j as i/j. We number the twelve versions as follows:             //
		//                                                                           //
		//           |   edge 0     edge 1     edge 2                                //
		//   --------|--------------------------------                               //
		//    face 0 |   0 (0/1)    4 (0/3)    8 (0/2)                               //
		//    face 1 |   1 (1/2)    5 (1/3)    9 (1/0)                               //
		//    face 2 |   2 (2/3)    6 (2/1)   10 (2/0)                               //
		//    face 3 |   3 (3/0)    7 (3/1)   11 (3/2)                               //
		//                                                                           //
		// Similarly, navigation and manipulation in a (boundary) triangulation are  //
		// done by using handles of triangles. Each handle is a pair (s, v), where s //
		// is a pointer to a triangle, and v is a version in the range from 0 to 5.  //
		// Each version corresponds to a directed edge of this triangle.             //
		//                                                                           //
		// Number the three vertices of a triangle from 0 to 2 (according to their   //
		// storage in the data structure). Give each edge the same index as the node //
		// opposite it in the triangle. The six versions of a triangle are:          //
		//                                                                           //
		//                 | edge 0   edge 1   edge 2                                //
		//  ---------------|--------------------------                               //
		//   ccw orieation |   0        2        4                                   //
		//    cw orieation |   1        3        5                                   //
		//                                                                           //
		// In the following, a 'triface' is a handle of tetrahedron, and a 'face' is //
		// a handle of a triangle.                                                   //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class triface
		{
		public:
			tetrahedron *tet;
			int ver; // Range from 0 to 11.
			triface() : tet(0), ver(0) {}
			triface& operator=(const triface& t) {
				tet = t.tet; ver = t.ver;
				return *this;
			}
		};

		class face
		{
		public:
			shellface *sh;
			int shver; // Range from 0 to 5.
			face() : sh(0), shver(0) {}
			face& operator=(const face& s) {
				sh = s.sh; shver = s.shver;
				return *this;
			}
		};

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Arraypool                                                                 //
		//                                                                           //
		// A dynamic linear array. (It is written by J. Shewchuk)                    //
		//                                                                           //
		// Each arraypool contains an array of pointers to a number of blocks.  Each //
		// block contains the same fixed number of objects.  Each index of the array //
		// addresses a particular object in the pool. The most significant bits add- //
		// ress the index of the block containing the object. The less significant   //
		// bits address this object within the block.                                //
		//                                                                           //
		// 'objectbytes' is the size of one object in blocks; 'log2objectsperblock'  //
		// is the base-2 logarithm of 'objectsperblock'; 'objects' counts the number //
		// of allocated objects; 'totalmemory' is the total memory in bytes.         //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class arraypool
		{
		public:

			int objectbytes;
			int objectsperblock;
			int log2objectsperblock;
			int objectsperblockmark;
			int toparraylen;
			char **toparray;
			long objects;
			unsigned long totalmemory;

			void restart();
			void poolinit(int sizeofobject, int log2objperblk);
			char* getblock(int objectindex);
			void* lookup(int objectindex);
			int newindex(void **newptr);

			arraypool(int sizeofobject, int log2objperblk);
			~arraypool();
		};

		// fastlookup() -- A fast, unsafe operation. Return the pointer to the object
		//   with a given index.  Note: The object's block must have been allocated,
		//   i.e., by the function newindex().

#define fastlookup(pool, index) \
	(void *)((pool)->toparray[(index) >> (pool)->log2objectsperblock] + \
	((index)& (pool)->objectsperblockmark) * (pool)->objectbytes)

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Memorypool                                                                //
		//                                                                           //
		// A structure for memory allocation. (It is written by J. Shewchuk)         //
		//                                                                           //
		// firstblock is the first block of items. nowblock is the block from which  //
		//   items are currently being allocated. nextitem points to the next slab   //
		//   of free memory for an item. deaditemstack is the head of a linked list  //
		//   (stack) of deallocated items that can be recycled.  unallocateditems is //
		//   the number of items that remain to be allocated from nowblock.          //
		//                                                                           //
		// Traversal is the process of walking through the entire list of items, and //
		//   is separate from allocation.  Note that a traversal will visit items on //
		//   the "deaditemstack" stack as well as live items.  pathblock points to   //
		//   the block currently being traversed.  pathitem points to the next item  //
		//   to be traversed.  pathitemsleft is the number of items that remain to   //
		//   be traversed in pathblock.                                              //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class memorypool
		{
		public:

			void **firstblock, **nowblock;
			void *nextitem;
			void *deaditemstack;
			void **pathblock;
			void *pathitem;
			int  alignbytes;
			int  itembytes, itemwords;
			int  itemsperblock;
			long items, maxitems;
			int  unallocateditems;
			int  pathitemsleft;

			memorypool();
			memorypool(int, int, int, int);
			~memorypool();

			void poolinit(int, int, int, int);
			void restart();
			void *alloc();
			void dealloc(void*);
			void traversalinit();
			void *traverse();
		};

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// badface                                                                   //
		//                                                                           //
		// Despite of its name, a 'badface' can be used to represent one of the      //
		// following objects:                                                        //
		//   - a face of a tetrahedron which is (possibly) non-Delaunay;             //
		//   - an encroached subsegment or subface;                                  //
		//   - a bad-quality tetrahedron, i.e, has too large radius-edge ratio;      //
		//   - a sliver, i.e., has good radius-edge ratio but nearly zero volume;    //
		//   - a recently flipped face (saved for undoing the flip later).           //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class badface {
		public:
			triface tt;
			face ss;
			REAL key, cent[6];  // circumcenter or cos(dihedral angles) at 6 edges.
			point forg, fdest, fapex, foppo, noppo;
			badface *nextitem;
			badface() : key(0), forg(0), fdest(0), fapex(0), foppo(0), noppo(0),
				nextitem(0) {}
		};

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// insertvertexflags                                                         //
		//                                                                           //
		// A collection of flags that pass to the routine insertvertex().            //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class insertvertexflags
		{

		public:

			int iloc;  // input/output.
			int bowywat, lawson;
			int splitbdflag, validflag, respectbdflag;
			int rejflag, chkencflag, cdtflag;
			int assignmeshsize;
			int sloc, sbowywat;

			// Used by Delaunay refinement.
			int refineflag; // 0, 1, 2, 3
			triface refinetet;
			face refinesh;
			int smlenflag; // for useinsertradius.
			REAL smlen; // for useinsertradius.
			point parentpt;

			insertvertexflags() {
				iloc = bowywat = lawson = 0;
				splitbdflag = validflag = respectbdflag = 0;
				rejflag = chkencflag = cdtflag = 0;
				assignmeshsize = 0;
				sloc = sbowywat = 0;

				refineflag = 0;
				refinetet.tet = NULL;
				refinesh.sh = NULL;
				smlenflag = 0;
				smlen = 0.0;
			}
		};

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// flipconstraints                                                           //
		//                                                                           //
		// A structure of a collection of data (options and parameters) which pass   //
		// to the edge flip function flipnm().                                       //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class flipconstraints {

		public:

			// Elementary flip flags.
			int enqflag; // (= flipflag)
			int chkencflag;

			// Control flags
			int unflip;  // Undo the performed flips.
			int collectnewtets; // Collect the new tets created by flips.
			int collectencsegflag;

			// Optimization flags.
			int remove_ndelaunay_edge; // Remove a non-Delaunay edge.
			REAL bak_tetprism_vol; // The value to be minimized.
			REAL tetprism_vol_sum;
			int remove_large_angle; // Remove a large dihedral angle at edge.
			REAL cosdihed_in; // The input cosine of the dihedral angle (> 0).
			REAL cosdihed_out; // The improved cosine of the dihedral angle.

			// Boundary recovery flags.
			int checkflipeligibility;
			point seg[2];  // A constraining edge to be recovered.
			point fac[3];  // A constraining face to be recovered.
			point remvert; // A vertex to be removed.


			flipconstraints() {
				enqflag = 0;
				chkencflag = 0;

				unflip = 0;
				collectnewtets = 0;
				collectencsegflag = 0;

				remove_ndelaunay_edge = 0;
				bak_tetprism_vol = 0.0;
				tetprism_vol_sum = 0.0;
				remove_large_angle = 0;
				cosdihed_in = 0.0;
				cosdihed_out = 0.0;

				checkflipeligibility = 0;
				seg[0] = NULL;
				fac[0] = NULL;
				remvert = NULL;
			}
		};

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// optparameters                                                             //
		//                                                                           //
		// Optimization options and parameters.                                      //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		class optparameters
		{
		public:
			// The one of goals of optimization.
			int max_min_volume;      // Maximize the minimum volume.
			int max_min_aspectratio; // Maximize the minimum aspect ratio.
			int min_max_dihedangle;  // Minimize the maximum dihedral angle. 

			// The initial and improved value.
			REAL initval, imprval;

			int numofsearchdirs;
			REAL searchstep;
			int maxiter;  // Maximum smoothing iterations (disabled by -1).
			int smthiter; // Performed iterations.


			optparameters() {
				max_min_volume = 0;
				max_min_aspectratio = 0;
				min_max_dihedangle = 0;

				initval = imprval = 0.0;

				numofsearchdirs = 10;
				searchstep = 0.01;
				maxiter = -1;   // Unlimited smoothing iterations.
				smthiter = 0;

			}
		};


		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Labels (enumeration declarations) used by TetGen.                         //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// Labels that signify the type of a vertex. 
		enum verttype
		{
			UNUSEDVERTEX, DUPLICATEDVERTEX, RIDGEVERTEX, ACUTEVERTEX,
			FACETVERTEX, VOLVERTEX, FREESEGVERTEX, FREEFACETVERTEX,
			FREEVOLVERTEX, NREGULARVERTEX, DEADVERTEX
		};

		// Labels that signify the result of triangle-triangle intersection test.
		enum interresult
		{
			DISJOINT, INTERSECT, SHAREVERT, SHAREEDGE, SHAREFACE,
			TOUCHEDGE, TOUCHFACE, ACROSSVERT, ACROSSEDGE, ACROSSFACE,
			COLLISIONFACE, ACROSSSEG, ACROSSSUB
		};

		// Labels that signify the result of point location.
		enum locateresult
		{
			UNKNOWN, OUTSIDE, INTETRAHEDRON, ONFACE, ONEDGE, ONVERTEX,
			ENCVERTEX, ENCSEGMENT, ENCSUBFACE, NEARVERTEX, NONREGULAR,
			INSTAR, BADELEMENT
		};

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Variables of TetGen                                                       //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// Pointer to the input data (a set of nodes, a PLC, or a mesh).
		tetgenio *in, *addin;

		// Pointer to the switches and parameters.
		tetgenbehavior b;

		// Pointer to a background mesh (contains size specification map).
		tetgenmesh *bgm;

		// Memorypools to store mesh elements (points, tetrahedra, subfaces, and
		//   segments) and extra pointers between tetrahedra, subfaces, and segments.
		memorypool *tetrahedrons, *subfaces, *subsegs, *points;
		memorypool *tet2subpool, *tet2segpool;

		// Memorypools to store bad-quality (or encroached) elements.
		memorypool *badtetrahedrons, *badsubfacs, *badsubsegs;

		// A memorypool to store faces to be flipped.
		memorypool *flippool;
		arraypool *unflipqueue;
		badface *flipstack;

		// Arrays used for point insertion (the Bowyer-Watson algorithm).
		arraypool *cavetetlist, *cavebdrylist, *caveoldtetlist;
		arraypool *cavetetshlist, *cavetetseglist, *cavetetvertlist;
		arraypool *caveencshlist, *caveencseglist;
		arraypool *caveshlist, *caveshbdlist, *cavesegshlist;

		// Stacks used for CDT construction and boundary recovery.
		arraypool *subsegstack, *subfacstack, *subvertstack;

		// Arrays of encroached segments and subfaces (for mesh refinement).
		arraypool *encseglist, *encshlist;

		// The map between facets to their vertices (for mesh refinement).
		int *idx2facetlist;
		point *facetverticeslist;

		// The map between segments to their endpoints (for mesh refinement).
		point *segmentendpointslist;

		// The infinite vertex.
		point dummypoint;
		// The recently visited tetrahedron, subface.
		triface recenttet;
		face recentsh;

		// PI is the ratio of a circle's circumference to its diameter.
		static const REAL PI;

		// Array (size = numberoftetrahedra * 6) for storing high-order nodes of
		//   tetrahedra (only used when -o2 switch is selected).
		point *highordertable;

		// Various variables.
		int numpointattrib;						// Number of point attributes.
		int numelemattrib;                     // Number of tetrahedron attributes.
		int sizeoftensor;                     // Number of REALs per metric tensor.
		int pointmtrindex;           // Index to find the metric tensor of a point.
		int pointparamindex;       // Index to find the u,v coordinates of a point.
		int point2simindex;         // Index to find a simplex adjacent to a point.
		int pointmarkindex;            // Index to find boundary marker of a point.
		int elemattribindex;          // Index to find attributes of a tetrahedron.
		int volumeboundindex;       // Index to find volume bound of a tetrahedron.
		int elemmarkerindex;              // Index to find marker of a tetrahedron.
		int shmarkindex;             // Index to find boundary marker of a subface.
		int areaboundindex;               // Index to find area bound of a subface.
		int checksubsegflag;   // Are there segments in the tetrahedralization yet?
		int checksubfaceflag;  // Are there subfaces in the tetrahedralization yet?
		int checkconstraints;  // Are there variant (node, seg, facet) constraints?
		int nonconvex;                               // Is current mesh non-convex?
		int autofliplinklevel;        // The increase of link levels, default is 1.
		int useinsertradius;       // Save the insertion radius for Steiner points.
		long samples;               // Number of random samples for point location.
		unsigned long randomseed;                    // Current random number seed.

		REAL cosmaxdihed, cosmindihed;    // The cosine values of max/min dihedral.
		REAL cossmtdihed;     // The cosine value of a bad dihedral to be smoothed.
		REAL cosslidihed;      // The cosine value of the max dihedral of a sliver.
		REAL minfaceang, minfacetdihed;     // The minimum input (dihedral) angles.
		REAL tetprism_vol_sum;   // The total volume of tetrahedral-prisms (in 4D).
		REAL longest;                          // The longest possible edge length.
		REAL xmax, xmin, ymax, ymin, zmax, zmin;         // Bounding box of points.

		// Counters.
		long insegments;                               // Number of input segments.
		long hullsize;                        // Number of exterior boundary faces.
		long meshedges;                                    // Number of mesh edges.
		long meshhulledges;                       // Number of boundary mesh edges.
		long steinerleft;                 // Number of Steiner points not yet used.
		long dupverts;                            // Are there duplicated vertices?
		long unuverts;                                // Are there unused vertices?
		long nonregularcount;                    // Are there non-regular vertices?
		long st_segref_count, st_facref_count, st_volref_count;  // Steiner points.
		long fillregioncount, cavitycount, cavityexpcount;
		long flip14count, flip26count, flipn2ncount;
		long flip23count, flip32count, flip44count, flip41count;
		long flip31count, flip22count;
		unsigned long totalworkmemory;      // Total memory used by working arrays.


		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh manipulation primitives                                              //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// Fast lookup tables for mesh manipulation primitives.
		static int bondtbl[12][12], fsymtbl[12][12];
		static int esymtbl[12], enexttbl[12], eprevtbl[12];
		static int enextesymtbl[12], eprevesymtbl[12];
		static int eorgoppotbl[12], edestoppotbl[12];
		static int facepivot1[12], facepivot2[12][12];
		static int orgpivot[12], destpivot[12], apexpivot[12], oppopivot[12];
		static int tsbondtbl[12][6], stbondtbl[12][6];
		static int tspivottbl[12][6], stpivottbl[12][6];
		static int ver2edge[12], edge2ver[6], epivot[12];
		static int sorgpivot[6], sdestpivot[6], sapexpivot[6];
		static int snextpivot[6];

		void inittables();

		// Primitives for tetrahedra.
		tetrahedron encode(triface& t);
		tetrahedron encode2(tetrahedron* ptr, int ver);
		void decode(tetrahedron ptr, triface& t);
		void bond(triface& t1, triface& t2);
		void dissolve(triface& t);
		void esym(triface& t1, triface& t2);
		void esymself(triface& t);
		void enext(triface& t1, triface& t2);
		void enextself(triface& t);
		void eprev(triface& t1, triface& t2);
		void eprevself(triface& t);
		void enextesym(triface& t1, triface& t2);
		void enextesymself(triface& t);
		void eprevesym(triface& t1, triface& t2);
		void eprevesymself(triface& t);
		void eorgoppo(triface& t1, triface& t2);
		void eorgoppoself(triface& t);
		void edestoppo(triface& t1, triface& t2);
		void edestoppoself(triface& t);
		void fsym(triface& t1, triface& t2);
		void fsymself(triface& t);
		void fnext(triface& t1, triface& t2);
		void fnextself(triface& t);
		point org(triface& t);
		point dest(triface& t);
		point apex(triface& t);
		point oppo(triface& t);
		void setorg(triface& t, point p);
		void setdest(triface& t, point p);
		void setapex(triface& t, point p);
		void setoppo(triface& t, point p);
		REAL elemattribute(tetrahedron* ptr, int attnum);
		void setelemattribute(tetrahedron* ptr, int attnum, REAL value);
		REAL volumebound(tetrahedron* ptr);
		void setvolumebound(tetrahedron* ptr, REAL value);
		int  elemindex(tetrahedron* ptr);
		void setelemindex(tetrahedron* ptr, int value);
		int  elemmarker(tetrahedron* ptr);
		void setelemmarker(tetrahedron* ptr, int value);
		void infect(triface& t);
		void uninfect(triface& t);
		bool infected(triface& t);
		void marktest(triface& t);
		void unmarktest(triface& t);
		bool marktested(triface& t);
		void markface(triface& t);
		void unmarkface(triface& t);
		bool facemarked(triface& t);
		void markedge(triface& t);
		void unmarkedge(triface& t);
		bool edgemarked(triface& t);
		void marktest2(triface& t);
		void unmarktest2(triface& t);
		bool marktest2ed(triface& t);
		int  elemcounter(triface& t);
		void setelemcounter(triface& t, int value);
		void increaseelemcounter(triface& t);
		void decreaseelemcounter(triface& t);
		bool ishulltet(triface& t);
		bool isdeadtet(triface& t);

		// Primitives for subfaces and subsegments.
		void sdecode(shellface sptr, face& s);
		shellface sencode(face& s);
		shellface sencode2(shellface *sh, int shver);
		void spivot(face& s1, face& s2);
		void spivotself(face& s);
		void sbond(face& s1, face& s2);
		void sbond1(face& s1, face& s2);
		void sdissolve(face& s);
		point sorg(face& s);
		point sdest(face& s);
		point sapex(face& s);
		void setsorg(face& s, point pointptr);
		void setsdest(face& s, point pointptr);
		void setsapex(face& s, point pointptr);
		void sesym(face& s1, face& s2);
		void sesymself(face& s);
		void senext(face& s1, face& s2);
		void senextself(face& s);
		void senext2(face& s1, face& s2);
		void senext2self(face& s);
		REAL areabound(face& s);
		void setareabound(face& s, REAL value);
		int shellmark(face& s);
		void setshellmark(face& s, int value);
		void sinfect(face& s);
		void suninfect(face& s);
		bool sinfected(face& s);
		void smarktest(face& s);
		void sunmarktest(face& s);
		bool smarktested(face& s);
		void smarktest2(face& s);
		void sunmarktest2(face& s);
		bool smarktest2ed(face& s);
		void smarktest3(face& s);
		void sunmarktest3(face& s);
		bool smarktest3ed(face& s);
		void setfacetindex(face& f, int value);
		int  getfacetindex(face& f);

		// Primitives for interacting tetrahedra and subfaces.
		void tsbond(triface& t, face& s);
		void tsdissolve(triface& t);
		void stdissolve(face& s);
		void tspivot(triface& t, face& s);
		void stpivot(face& s, triface& t);

		// Primitives for interacting tetrahedra and segments.
		void tssbond1(triface& t, face& seg);
		void sstbond1(face& s, triface& t);
		void tssdissolve1(triface& t);
		void sstdissolve1(face& s);
		void tsspivot1(triface& t, face& s);
		void sstpivot1(face& s, triface& t);

		// Primitives for interacting subfaces and segments.
		void ssbond(face& s, face& edge);
		void ssbond1(face& s, face& edge);
		void ssdissolve(face& s);
		void sspivot(face& s, face& edge);

		// Primitives for points.
		int  pointmark(point pt);
		void setpointmark(point pt, int value);
		enum verttype pointtype(point pt);
		void setpointtype(point pt, enum verttype value);
		int  pointgeomtag(point pt);
		void setpointgeomtag(point pt, int value);
		REAL pointgeomuv(point pt, int i);
		void setpointgeomuv(point pt, int i, REAL value);
		void pinfect(point pt);
		void puninfect(point pt);
		bool pinfected(point pt);
		void pmarktest(point pt);
		void punmarktest(point pt);
		bool pmarktested(point pt);
		void pmarktest2(point pt);
		void punmarktest2(point pt);
		bool pmarktest2ed(point pt);
		void pmarktest3(point pt);
		void punmarktest3(point pt);
		bool pmarktest3ed(point pt);
		tetrahedron point2tet(point pt);
		void setpoint2tet(point pt, tetrahedron value);
		shellface point2sh(point pt);
		void setpoint2sh(point pt, shellface value);
		point point2ppt(point pt);
		void setpoint2ppt(point pt, point value);
		tetrahedron point2bgmtet(point pt);
		void setpoint2bgmtet(point pt, tetrahedron value);
		void setpointinsradius(point pt, REAL value);
		REAL getpointinsradius(point pt);

		// Advanced primitives.
		void point2tetorg(point pt, triface& t);
		void point2shorg(point pa, face& s);
		point farsorg(face& seg);
		point farsdest(face& seg);

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		//  Memory managment                                                         //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		void tetrahedrondealloc(tetrahedron*);
		tetrahedron *tetrahedrontraverse();
		tetrahedron *alltetrahedrontraverse();
		void shellfacedealloc(memorypool*, shellface*);
		shellface *shellfacetraverse(memorypool*);
		void pointdealloc(point);
		point pointtraverse();

		void makeindex2pointmap(point*&);
		void makepoint2submap(memorypool*, int*&, face*&);
		void maketetrahedron(triface*);
		void makeshellface(memorypool*, face*);
		void makepoint(point*, enum verttype);

		void initializepools();

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Advanced geometric predicates and calculations                            //
		//                                                                           //
		// TetGen uses a simplified symbolic perturbation scheme from Edelsbrunner,  //
		// et al [*].  Hence the point-in-sphere test never returns a zero. The idea //
		// is to perturb the weights of vertices in the fourth dimension.  TetGen    //
		// uses the indices of the vertices decide the amount of perturbation. It is //
		// implemented in the routine insphere_s().
		//                                                                           //
		// The routine tri_edge_test() determines whether or not a triangle and an   //
		// edge intersect in 3D. If they intersect, their intersection type is also  //
		// reported. This test is a combination of n 3D orientation tests (n is bet- //
		// ween 3 and 9). It uses the robust orient3d() test to make the branch dec- //
		// isions.  The routine tri_tri_test() determines whether or not two triang- //
		// les intersect in 3D. It also uses the robust orient3d() test.             //
		//                                                                           //
		// There are a number of routines to calculate geometrical quantities, e.g., //
		// circumcenters, angles, dihedral angles, face normals, face areas, etc.    //
		// They are so far done by the default floating-point arithmetics which are  //
		// non-robust. They should be improved in the future.                        //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// Symbolic perturbations (robust)
		REAL insphere_s(REAL*, REAL*, REAL*, REAL*, REAL*);
		REAL orient4d_s(REAL*, REAL*, REAL*, REAL*, REAL*,
			REAL, REAL, REAL, REAL, REAL);

		// Triangle-edge intersection test (robust)
		int tri_edge_2d(point, point, point, point, point, point, int, int*, int*);
		int tri_edge_tail(point, point, point, point, point, point, REAL, REAL, int,
			int*, int*);
		int tri_edge_test(point, point, point, point, point, point, int, int*, int*);

		// Triangle-triangle intersection test (robust)
		int tri_edge_inter_tail(point, point, point, point, point, REAL, REAL);
		int tri_tri_inter(point, point, point, point, point, point);

		// Linear algebra functions
		REAL dot(REAL* v1, REAL* v2);
		void cross(REAL* v1, REAL* v2, REAL* n);
		bool lu_decmp(REAL lu[4][4], int n, int* ps, REAL* d, int N);
		void lu_solve(REAL lu[4][4], int n, int* ps, REAL* b, int N);

		// An embedded 2-dimensional geometric predicate (non-robust)
		REAL incircle3d(point pa, point pb, point pc, point pd);

		// Geometric calculations (non-robust)
		REAL orient3dfast(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
		REAL norm2(REAL x, REAL y, REAL z);
		REAL distance(REAL* p1, REAL* p2);
		void facenormal(point pa, point pb, point pc, REAL *n, int pivot, REAL *lav);
		REAL shortdistance(REAL* p, REAL* e1, REAL* e2);
		REAL triarea(REAL* pa, REAL* pb, REAL* pc);
		REAL interiorangle(REAL* o, REAL* p1, REAL* p2, REAL* n);
		void projpt2edge(REAL* p, REAL* e1, REAL* e2, REAL* prj);
		void projpt2face(REAL* p, REAL* f1, REAL* f2, REAL* f3, REAL* prj);
		REAL facedihedral(REAL* pa, REAL* pb, REAL* pc1, REAL* pc2);
		bool tetalldihedral(point, point, point, point, REAL*, REAL*, REAL*);
		void tetallnormal(point, point, point, point, REAL N[4][3], REAL* volume);
		REAL tetaspectratio(point, point, point, point);
		bool circumsphere(REAL*, REAL*, REAL*, REAL*, REAL* cent, REAL* radius);
		bool orthosphere(REAL*, REAL*, REAL*, REAL*, REAL, REAL, REAL, REAL, REAL*, REAL*);
		void planelineint(REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*);
		int linelineint(REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*, REAL*);
		REAL tetprismvol(REAL* pa, REAL* pb, REAL* pc, REAL* pd);
		bool calculateabovepoint(arraypool*, point*, point*, point*);
		void calculateabovepoint4(point, point, point, point);

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Local mesh transformations                                                //
		//                                                                           //
		// A local transformation replaces a small set of tetrahedra with another    //
		// set of tetrahedra which fills the same space and the same boundaries.     //
		//   In 3D, the most simplest local transformations are the elementary flips //
		// performed within the convex hull of five vertices: 2-to-3, 3-to-2, 1-to-4,//
		// and 4-to-1 flips,  where the numbers indicate the number of tetrahedra    //
		// before and after each flip.  The 1-to-4 and 4-to-1 flip involve inserting //
		// or deleting a vertex, respectively.                                       //
		//   There are complex local transformations which can be decomposed as a    //
		// combination of elementary flips. For example,a 4-to-4 flip which replaces //
		// two coplanar edges can be regarded by a 2-to-3 flip and a 3-to-2 flip.    //
		// Note that the first 2-to-3 flip will temporarily create a degenerate tet- //
		// rahedron which is removed immediately by the followed 3-to-2 flip.  More  //
		// generally, a n-to-m flip, where n > 3, m = (n - 2) * 2, which removes an  //
		// edge can be done by first performing a sequence of (n - 3) 2-to-3 flips   //
		// followed by a 3-to-2 flip.                                                //
		//                                                                           //
		// The routines flip23(), flip32(), and flip41() perform the three element-  //
		// ray flips. The flip14() is available inside the routine insertpoint().    //
		//                                                                           //
		// The routines flipnm() and flipnm_post() implement a generalized edge flip //
		// algorithm which uses a combination of elementary flips.                   //
		//                                                                           //
		// The routine insertpoint() implements a variant of Bowyer-Watson's cavity  //
		// algorithm to insert a vertex. It works for arbitrary tetrahedralization,  //
		// either Delaunay, or constrained Delaunay, or non-Delaunay.                //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// The elementary flips.
		void flip23(triface*, int, flipconstraints* fc);
		void flip32(triface*, int, flipconstraints* fc);
		void flip41(triface*, int, flipconstraints* fc);

		// A generalized edge flip.
		int flipnm(triface*, int n, int level, int, flipconstraints* fc);
		int flipnm_post(triface*, int n, int nn, int, flipconstraints* fc);

		// Point insertion.
		int  insertpoint(point, triface*, face*, face*, insertvertexflags*);
		void insertpoint_abort(face*, insertvertexflags*);

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Delaunay tetrahedralization                                               //
		//                                                                           //
		// The routine incrementaldelaunay() implemented two incremental algorithms  //
		// for constructing Delaunay tetrahedralizations (DTs):  the Bowyer-Watson   //
		// (B-W) algorithm and the incremental flip algorithm of Edelsbrunner and    //
		// Shah, "Incremental topological flipping works for regular triangulation," //
		// Algorithmica, 15:233-241, 1996.                                           //
		//                                                                           //
		// The routine incrementalflip() implements the flip algorithm of [Edelsbru- //
		// nner and Shah, 1996].  It flips a queue of locally non-Delaunay faces (in //
		// an arbitrary order).  The success is guaranteed when the Delaunay tetrah- //
		// edralization is constructed incrementally by adding one vertex at a time. //
		//                                                                           //
		// The routine locate() finds a tetrahedron contains a new point in current  //
		// DT.  It uses a simple stochastic walk algorithm: starting from an arbitr- //
		// ary tetrahedron in DT, it finds the destination by visit one tetrahedron  //
		// at a time, randomly chooses a tetrahedron if there are more than one      //
		// choices. This algorithm terminates due to Edelsbrunner's acyclic theorem. //
		//   Choose a good starting tetrahedron is crucial to the speed of the walk. //
		// TetGen originally uses the "jump-and-walk" algorithm of Muecke, E.P.,     //
		// Saias, I., and Zhu, B. "Fast Randomized Point Location Without Preproces- //
		// sing." In Proceedings of the 12th ACM Symposium on Computational Geometry,//
		// 274-283, 1996.  It first randomly samples several tetrahedra in the DT    //
		// and then choosing the closet one to start walking.                        //
		//   The above algorithm slows download dramatically as the number of points //
		// grows -- reported in Amenta, N., Choi, S. and Rote, G., "Incremental      //
		// construction con {BRIO}," In Proceedings of 19th ACM Symposium on         //
		// Computational Geometry, 211-219, 2003.  On the other hand, Liu and        //
		// Snoeyink showed that the point location can be made in constant time if   //
		// the points are pre-sorted so that the nearby points in space have nearby  //
		// indices, then adding the points in this order. They sorted the points     //
		// along the 3D Hilbert curve.                                               //
		//                                                                           //
		// The routine hilbert_sort3() sorts a set of 3D points along the 3D Hilbert //
		// curve. It recursively splits a point set according to the Hilbert indices //
		// mapped to the subboxes of the bounding box of the point set.              //
		//   The Hilbert indices is calculated by Butz's algorithm in 1971.  A nice  //
		// exposition of this algorithm can be found in the paper of Hamilton, C.,   //
		// "Compact Hilbert Indices", Technical Report CS-2006-07, Computer Science, //
		// Dalhousie University, 2006 (the Section 2). My implementation also refer- //
		// enced Steven Witham's implementation of "Hilbert walk" (hopefully, it is  //
		// still available at: http://www.tiac.net/~sw/2008/10/Hilbert/).            //
		//                                                                           //
		// TetGen sorts the points using the method in the paper of Boissonnat,J.-D.,//
		// Devillers, O. and Hornus, S. "Incremental Construction of the Delaunay    //
		// Triangulation and the Delaunay Graph in Medium Dimension," In Proceedings //
		// of the 25th ACM Symposium on Computational Geometry, 2009.                //
		//   It first randomly sorts the points into subgroups using the Biased Rand-//
		// omized Insertion Ordering (BRIO) of Amenta et al 2003, then sorts the     //
		// points in each subgroup along the 3D Hilbert curve.  Inserting points in  //
		// this order ensures a randomized "sprinkling" of the points over the       //
		// domain, while sorting of each subset ensures locality.                    //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		void transfernodes();

		// Point sorting.
		int  transgc[8][3][8], tsb1mod3[8];
		void hilbert_init(int n);
		int  hilbert_split(point* vertexarray, int arraysize, int gc0, int gc1,
			REAL, REAL, REAL, REAL, REAL, REAL);
		void hilbert_sort3(point* vertexarray, int arraysize, int e, int d,
			REAL, REAL, REAL, REAL, REAL, REAL, int depth);
		void brio_multiscale_sort(point*, int, int threshold, REAL ratio, int* depth);

		// Point location.
		unsigned long randomnation(unsigned int choices);
		void randomsample(point searchpt, triface *searchtet);
		enum locateresult locate(point searchpt, triface *searchtet);

		// Incremental flips.
		void flippush(badface*&, triface*);
		int  incrementalflip(point newpt, int, flipconstraints *fc);

		// Incremental Delaunay construction.
		void initialdelaunay(point pa, point pb, point pc, point pd);
		void incrementaldelaunay(clock_t&);

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Surface triangulation                                                     //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		void flipshpush(face*);
		void flip22(face*, int, int);
		void flip31(face*, int);
		long lawsonflip();
		int sinsertvertex(point newpt, face*, face*, int iloc, int bowywat, int);
		int sremovevertex(point delpt, face*, face*, int lawson);

		enum locateresult slocate(point, face*, int, int, int);
		enum interresult sscoutsegment(face*, point);
		void scarveholes(int, REAL*);
		void triangulate(int, arraypool*, arraypool*, int, REAL*);

		void unifysubfaces(face*, face*);
		void unifysegments();
		void mergefacets();
		void identifypscedges(point*);
		void meshsurface();

		void interecursive(shellface** subfacearray, int arraysize, int axis,
			REAL, REAL, REAL, REAL, REAL, REAL, int* internum);
		void detectinterfaces();

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Constrained Delaunay tetrahedralization                                   //
		//                                                                           //
		// A constrained Delaunay tetrahedralization (CDT) is a variation of a Dela- //
		// unay tetrahedralization (DT) that is constrained to respect the boundary  //
		// of a 3D PLC (domain). In a CDT of a 3D PLC, every vertex or edge of the   //
		// PLC is also a vertex or an edge of the CDT, every polygon of the PLC is a //
		// union of triangles of the CDT. A crucial difference between a CDT and a   //
		// DT is that triangles in the PLC's polygons are not required to be locally //
		// Delaunay, which frees the CDT to better respect the PLC's polygons. CDTs  //
		// have optimal properties similar to those of DTs.                          //
		//                                                                           //
		// Steiner Points and Steiner CDTs. It is known that even a simple 3D polyh- //
		// edron may not have a tetrahedralization which only uses its own vertices. //
		// Some extra points, so-called "Steiner points" are needed in order to form //
		// a tetrahedralization of such polyhedron.  It is true for tetrahedralizing //
		// a 3D PLC as well. A Steiner CDT of a 3D PLC is a CDT containing Steiner   //
		// points. The CDT algorithms of TetGen in general create Steiner CDTs.      //
		// Almost all of the Steiner points are added in the edges of the PLC. They  //
		// guarantee the existence of a CDT of the modified PLC.                     //
		//                                                                           //
		// The routine constraineddelaunay() starts from a DT of the vertices of a   //
		// PLC and creates a (Steiner) CDT of the PLC (including Steiner points). It //
		// is constructed by two steps, (1) segment recovery and (2) facet (polygon) //
		// recovery. Each step is accomplished by its own algorithm.                 //
		//                                                                           //
		// The routine delaunizesegments() implements the segment recovery algorithm //
		// of Si, H. and Gaertner, K. "Meshing Piecewise Linear Complexes by Constr- //
		// ained Delaunay Tetrahedralizations," In Proceedings of the 14th Internat- //
		// ional Meshing Roundtable, 147--163, 2005.  It adds Steiner points into    //
		// non-Delaunay segments until all subsegments appear together in a DT. The  //
		// running time of this algorithm is proportional to the number of added     //
		// Steiner points.                                                           //
		//                                                                           //
		// There are two incremental facet recovery algorithms: the cavity re-trian- //
		// gulation algorithm of Si, H. and Gaertner, K. "3D Boundary Recovery by    //
		// Constrained Delaunay Tetrahedralization," International Journal for Numer-//
		// ical Methods in Engineering, 85:1341-1364, 2011, and the flip algorithm   //
		// of Shewchuk, J. "Updating and Constructing Constrained Delaunay and       //
		// Constrained Regular Triangulations by Flips." In Proceedings of the 19th  //
		// ACM Symposium on Computational Geometry, 86-95, 2003.                     //
		//                                                                           //
		// It is guaranteed in theory, no Steiner point is needed in both algorithms //
		// However, a facet with non-coplanar vertices might cause the  additions of //
		// Steiner points. It is discussed in the paper of Si, H., and  Shewchuk, J.,//
		// "Incrementally Constructing and Updating Constrained Delaunay             //
		// Tetrahedralizations with Finite Precision Coordinates." In Proceedings of //
		// the 21th International Meshing Roundtable, 2012.                          //
		//                                                                           //
		// Our implementation of the facet recovery algorithms recover a "missing    //
		// region" at a time. Each missing region is a subset of connected interiors //
		// of a polygon. The routine formcavity() creates the cavity of crossing     //
		// tetrahedra of the missing region.                                         //
		//                                                                           //
		// The cavity re-triangulation algorithm is implemented by three subroutines,//
		// delaunizecavity(), fillcavity(), and carvecavity(). Since it may fail due //
		// to non-coplanar vertices, the subroutine restorecavity() is used to rest- //
		// ore the original cavity.                                                  //
		//                                                                           //
		// The routine flipinsertfacet() implements the flip algorithm. The subrout- //
		// ine flipcertify() is used to maintain the priority queue of flips.        // 
		//                                                                           //
		// The routine refineregion() is called when the facet recovery algorithm    //
		// fail to recover a missing region. It inserts Steiner points to refine the //
		// missing region. In order to avoid inserting Steiner points very close to  //
		// existing segments.  The classical encroachment rules of the Delaunay      //
		// refinement algorithm are used to choose the Steiner points.               //
		//                                                                           //
		// The routine constrainedfacets() does the facet recovery by using either   //
		// the cavity re-triangulation algorithm (default) or the flip algorithm. It //
		// results a CDT of the (modified) PLC (including Steiner points).           //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		void makesegmentendpointsmap();

		enum interresult finddirection(triface* searchtet, point endpt);
		enum interresult scoutsegment(point, point, triface*, point*, arraypool*);
		int  getsteinerptonsegment(face* seg, point refpt, point steinpt);
		void delaunizesegments();

		enum interresult scoutsubface(face* searchsh, triface* searchtet);
		void formregion(face*, arraypool*, arraypool*, arraypool*);
		int  scoutcrossedge(triface& crosstet, arraypool*, arraypool*);
		bool formcavity(triface*, arraypool*, arraypool*, arraypool*, arraypool*,
			arraypool*, arraypool*);

		// Facet recovery by cavity re-triangulation [Si and Gaertner 2011].
		void delaunizecavity(arraypool*, arraypool*, arraypool*, arraypool*,
			arraypool*, arraypool*);
		bool fillcavity(arraypool*, arraypool*, arraypool*, arraypool*,
			arraypool*, arraypool*, triface* crossedge);
		void carvecavity(arraypool*, arraypool*, arraypool*);
		void restorecavity(arraypool*, arraypool*, arraypool*, arraypool*);

		// Facet recovery by flips [Shewchuk 2003].
		void flipcertify(triface *chkface, badface **pqueue, point, point, point);
		void flipinsertfacet(arraypool*, arraypool*, arraypool*, arraypool*);

		bool fillregion(arraypool* missingshs, arraypool*, arraypool* newshs);

		int  insertpoint_cdt(point, triface*, face*, face*, insertvertexflags*,
			arraypool*, arraypool*, arraypool*, arraypool*,
			arraypool*, arraypool*);
		void refineregion(face&, arraypool*, arraypool*, arraypool*, arraypool*,
			arraypool*, arraypool*);

		void constrainedfacets();

		void constraineddelaunay(clock_t&);

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Constrained tetrahedralizations.                                          //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		int checkflipeligibility(int fliptype, point, point, point, point, point,
			int level, int edgepivot, flipconstraints* fc);

		int removeedgebyflips(triface*, flipconstraints*);
		int removefacebyflips(triface*, flipconstraints*);

		int recoveredgebyflips(point, point, triface*, int fullsearch);
		int add_steinerpt_in_schoenhardtpoly(triface*, int, int chkencflag);
		int add_steinerpt_in_segment(face*, int searchlevel);
		int addsteiner4recoversegment(face*, int);
		int recoversegments(arraypool*, int fullsearch, int steinerflag);

		int recoverfacebyflips(point, point, point, face*, triface*);
		int recoversubfaces(arraypool*, int steinerflag);

		int getvertexstar(int, point searchpt, arraypool*, arraypool*, arraypool*);
		int getedge(point, point, triface*);
		int reduceedgesatvertex(point startpt, arraypool* endptlist);
		int removevertexbyflips(point steinerpt);

		int suppressbdrysteinerpoint(point steinerpt);
		int suppresssteinerpoints();

		void recoverboundary(clock_t&);

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh reconstruction                                                       //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////
		void carveholes();

		void reconstructmesh();

		int  scoutpoint(point, triface*, int randflag);
		REAL getpointmeshsize(point, triface*, int iloc);
		void interpolatemeshsize();

		void insertconstrainedpoints(point *insertarray, int arylen, int rejflag);
		void insertconstrainedpoints(tetgenio *addio);

		void collectremovepoints(arraypool *remptlist);
		void meshcoarsening();

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh refinement                                                           //
		//                                                                           //
		// The purpose of mesh refinement is to obtain a tetrahedral mesh with well- //
		// -shaped tetrahedra and appropriate mesh size.  It is necessary to insert  //
		// new Steiner points to achieve this property. The questions are (1) how to //
		// choose the Steiner points? and (2) how to insert them?                    //
		//                                                                           //
		// Delaunay refinement is a technique first developed by Chew [1989] and     //
		// Ruppert [1993, 1995] to generate quality triangular meshes in the plane.  //
		// It provides guarantee on the smallest angle of the triangles.  Rupper's   //
		// algorithm guarantees that the mesh is size-optimal (to within a constant  //
		// factor) among all meshes with the same quality.                           //
		//   Shewchuk generalized Ruppert's algorithm into 3D in his PhD thesis      //
		// [Shewchuk 1997]. A short version of his algorithm appears in "Tetrahedral //
		// Mesh Generation by Delaunay Refinement," In Proceedings of the 14th ACM   //
		// Symposium on Computational Geometry, 86-95, 1998.  It guarantees that all //
		// tetrahedra of the output mesh have a "radius-edge ratio" (equivalent to   //
		// the minimal face angle) bounded. However, it does not remove slivers, a   //
		// type of very flat tetrahedra which can have no small face angles but have //
		// very small (and large) dihedral angles. Moreover, it may not terminate if //
		// the input PLC contains "sharp features", e.g., two edges (or two facets)  //
		// meet at an acute angle (or dihedral angle).                               //
		//                                                                           //
		// TetGen uses the basic Delaunay refinement scheme to insert Steiner points.//
		// While it always maintains a constrained Delaunay mesh.  The algorithm is  //
		// described in Si, H., "Adaptive Constrained Delaunay Mesh Generation,"     //
		// International Journal for Numerical Methods in Engineering, 75:856-880.   //
		// This algorithm always terminates and sharp features are easily preserved. //
		// The mesh has good quality (same as Shewchuk's Delaunay refinement algori- //
		// thm) in the bulk of the mesh domain. Moreover, it supports the generation //
		// of adaptive mesh according to a (isotropic) mesh sizing function.         //   
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		void makefacetverticesmap();
		int segsegadjacent(face *, face *);
		int segfacetadjacent(face *checkseg, face *checksh);
		int facetfacetadjacent(face *, face *);

		int checkseg4encroach(point pa, point pb, point checkpt);
		int checkseg4split(face *chkseg, point&, int&);
		int splitsegment(face *splitseg, point encpt, REAL, point, point, int, int);
		void repairencsegs(int chkencflag);

		void enqueuesubface(memorypool*, face*);
		int checkfac4encroach(point, point, point, point checkpt, REAL*, REAL*);
		int checkfac4split(face *chkfac, point& encpt, int& qflag, REAL *ccent);
		int splitsubface(face *splitfac, point, point, int qflag, REAL *ccent, int);
		void repairencfacs(int chkencflag);

		void enqueuetetrahedron(triface*);
		int checktet4split(triface *chktet, int& qflag, REAL *ccent);
		int splittetrahedron(triface* splittet, int qflag, REAL *ccent, int);
		void repairbadtets(int chkencflag);

		void delaunayrefinement();

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh optimization                                                         //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		long lawsonflip3d(flipconstraints *fc);
		void recoverdelaunay();

		int  gettetrahedron(point, point, point, point, triface *);
		long improvequalitybyflips();

		int  smoothpoint(point smtpt, arraypool*, int ccw, optparameters *opm);
		long improvequalitybysmoothing(optparameters *opm);

		int  splitsliver(triface *, REAL, int);
		long removeslivers(int);

		void optimizemesh();

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh check and statistics                                                 //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		// Mesh validations.
		int checkmesh(int topoflag);
		int checkshells();
		int checksegments();
		int checkdelaunay();
		int checkregular(int);
		int checkconforming(int);

		//  Mesh statistics.
		void printfcomma(unsigned long n);
		void qualitystatistics();
		void memorystatistics();
		void statistics();

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// Mesh output                                                               //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////

		void jettisonnodes();
		void highorder();
		void numberedges();
		void outnodes(tetgenio*);
		void outmetrics(tetgenio*);
		void outelements(tetgenio*);
		void outfaces(tetgenio*);
		void outhullfaces(tetgenio*);
		void outsubfaces(tetgenio*);
		void outedges(tetgenio*);
		void outsubsegments(tetgenio*);
		void outneighbors(tetgenio*);
		void outvoronoi(tetgenio*);
		void outsmesh(char*);
		void outmesh2medit(char*);
		void outmesh2vtk(char*);

		void outmesh2tecplot(char*);  //by gzz
		void outmesh2guo(const int& cpart);					//by gzz

		///////////////////////////////////////////////////////////////////////////////
		//                                                                           //
		// terminatetetgen()    Terminate TetGen with a given exit code.             //
		//                                                                           //
		///////////////////////////////////////////////////////////////////////////////
		static void terminatetetgen(int x);		//终止tetgen

	protected:
		void setvertices(triface& t, REAL* torg, REAL* tdest, REAL* tapex, REAL* toppo);
		// Quickly check if the handle (t, v) is a subface.
		bool issubface(const triface& t);
		void setshvertices(face& s, REAL* pa, REAL* pb, REAL* pc);

		// Quickly check if this subface is attached to a tetrahedron.
		bool isshtet(const face& s);

		// Quickly check if the edge is a subsegment.
		bool isshsubseg(const face& s);

		// Quickly check whether 't' is a segment or not.
		bool issubseg(const triface& t);
	};	

	/*//////////////////////////////////////////////////////////////////////////////
	//                                                                           //
	// Robust Geometric predicates                                               //
	//                                                                           //
	// Geometric predicates are simple tests of spatial relations of a set of d- //
	// dimensional points, such as the orientation test and the point-in-sphere  //
	// test. Each of these tests is performed by evaluating the sign of a deter- //
	// minant of a matrix whose entries are the coordinates of these points.  If //
	// the computation is performed by using the floating-point numbers, e.g.,   //
	// the single or double precision numbers in C/C++, roundoff error may cause //
	// an incorrect result. This may either lead to a wrong result or eventually //
	// lead to a failure of the program.  Computing the predicates exactly will  //
	// avoid the error and make the program robust.                              //
	//                                                                           //
	// The following routines are the robust geometric predicates for 3D orient- //
	// ation test and point-in-sphere test.  They were implemented by Shewchuk.  //
	// The source code are generously provided by him in the public domain,      //
	// http://www.cs.cmu.edu/~quake/robust.html. predicates.cxx is a C++ version //
	// of the original C code.                                                   //
	//                                                                           //
	// The original predicates of Shewchuk only use "dynamic filters", i.e., it  //
	// computes the error at run time step by step. TetGen first adds a "static  //
	// filter" in each predicate. It estimates the maximal possible error in all //
	// cases.  So it can safely and quickly answer many easy cases.              //
	//                                                                           //
	//////////////////////////////////////////////////////////////////////////////*/
	void exactinit(int, int, int, REAL, REAL, REAL);
	REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
	REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe);
	REAL orient4d(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe, REAL ah, REAL bh, REAL ch, REAL dh, REAL eh);
}